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Abstract 

We study a three-mode Hamiltonian modelling a heteronuclear molecular Bose- 
Einstein condensate. Two modes are associated with two distinguishable atomic 
constituents, which can combine to form a molecule represented by the third mode. 
Beginning with a semi-classical analogue of the model, we conduct an analysis to 
determine the phase space fixed points of the system. Bifurcations of the fixed 
points naturally separate the coupling parameter space into different regions. Two 
distinct scenarios are found, dependent on whether the imbalance between the num- 
ber operators for the atomic modes is zero or non-zero. This result suggests the 
ground-state properties of the model exhibit an unusual sensitivity on the atomic 
imbalance. We then test this finding for the quantum mechanical model. Specifi- 
cally we use Bethe ansatz methods, ground-state expectation values, the character 
of the quantum dynamics, and ground-state wavefunction overlaps to clarify the na- 
ture of the ground-state phases. The character of the transition is smoothed due to 
quantum fluctuations, but we may nonetheless identify the emergence of a quantum 
phase boundary in the limit of zero atomic imbalance. 

PACS: 02.30.1k, 03.65. Sq, 03.75.Nt 



1 Introduction 



The achievement of producing Bose-Einstein condenstates with ultracold dilute gases of 
atoms has seen a wealth of theoretical and experimental activity. One enticing prospect 
of Bose-Einstein condensates is that they may allow for a better understanding of the 
interface between classical and quantum mechnics, through the possibility of macroscopic 
Schrodinger cat states [1] and macroscopic quantum tunneling [2] . Another intriguing field 
of study is the chemistry of Bose-Einstein condensates, where the atomic constituents may 
form molecules through Feshbach resonances [3] or photoassociation [4]. A novel feature 
of a molecular Bose-Einstein condensate is that the atomic and molecular states can 
exist as a superposition [5], providing a chemical analogue of a Schrodinger cat state. In 
cases where the molecules are heteronuclear, the presence of a permanent electric dipole 
moment also opens the possibility for manipulating the condensate through electrostatic 
forces [6]. 

Since systems of Bose-Einstein condensates exist at ultracold temperatures, it is to 
be expected that significant insights into their behaviour can be obtained from studying 
their ground-state properties. From a general theoretical perspective there has been 
substantial progress in the understanding of quantum (i.e. ground-state) phases in many- 
body quantum systems, due largely to a cross fertilisation of ideas between the condensed 
matter theory and the quantum information theory communities. Much of this study has 
explored the relationship between entanglement and quantum criticality [7-9]. However 
other characterisations of quantum criticality have been sought too [10,11]. Recently the 
notion of wavefunction overlaps (also known as the fidelity), which is again common in 
quantum information theory, has been applied to the study of quantum phase transitions 
[12, 13]. An advantage of this approach is its universality, as it can be applied to any 
system independent of the choice of decomposition into subsystems. 

With the above points in mind here we analyse a simple, yet non-trivial, three-mode 
model describing a heteronuclear molecular condensate. Two modes are associated with 
two distinguishable atomic constituents, which can combine to form a molecule repre- 
sented by the third mode. Besides the interaction describing the interconversion of atoms 
and molecules, the Hamiltonian contains terms which are linear in the mode number op- 
erators (corresponding to external fields) and terms which are second-order in the mode 
number operators (corresponding to scattering interactions between atoms and molecules). 
We mainly concern ourselves with the ground-state properties of the model, with the aim 
of identifying the ground-state phases. We avoid taking the thermodynamic limit and 
restrict our analysis to finite systems, for reasons which will be discussed later. This in 
turn presents challenges in rigourously identifying quantum phases, since for finite sys- 
tems there are no singularities in physical quantities such as the ground state energy 
and its derivatives. However several recent works have addressed the issues of quantum 
phases in finite systems [14-17]. We mention that traditional techniques of renormali- 
sation group methods are not applicable to the model under consideration, due to the 
low number of degrees of freedom. Neither is the concept of symmetry breaking, as the 
model does not admit global symmetries, nor long-range order, as the model is in essence 
zero-dimensional. 

We start our analysis with a semi-classical treatment, following the approach of [18]. 
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Since the model with which we are dealing is integrable, the semi-classical many-body 
system can be reduced to a problem with a single degree of freedom. We study the 
phase space of this system, in particular determining the fixed points. It is found that for 
certain coupling parameters bifurcations of the fixed points occur, and we can determine a 
parameter space diagram which classifies the fixed points. An unexpected result is that the 
boundaries between the regions in parameter space are extremely sensitive on whether the 
number of constituent atoms is equal or not. Specifically, when the number of constituent 
atoms is equal (i.e. the atomic imbalance is zero) there is a spontaneous appearance 
of additional boundaries in the parameter space, some of which can be identified with 
bifurcations of the global minimum of the classical Hamiltonian. 

We next investigate the extent to which the classical behaviour influences the ground- 
state properties of the quantum system. Our first goal in the full quantum analysis is to 
derive an exact Bethe ansatz solution for the model. We use the Bethe ansatz solution 
to map the spectrum of the Hamiltonian into that of a one-body Schrodinger equation 
in one-dimensional. An advantage of this method is that it allows for an analysis of the 
finite system, following the ideas of [15], as the mapping to the one-body Schrodinger 
equation is not dependent on taking the thermodynamic limit of the original many-body 
system. The results of the analysis of the associated Schrodinger equation are in general 
agreement with the results obtained from the semi-classical treatment, supporting the 
picture of an additional phase boundary when the atomic imbalance is zero. However, 
due to quantum fluctuations, the emergence of the phase boundary is smooth rather than 
spontaneous. This property is apparent from a study of ground-state expectation values 
and quantum dynamics. 

In order to simply characterise the ground-state phases for the finite system, we finally 
define the notion of a quantum phase pre-transition in terms of wavefunction overlaps. 
Specifically, a quantum phase pre-transition is identified with each coupling for which the 
incremental ground-state wavefunction overlap is a local minimum. We numerically cal- 
culate these for several cases and discuss these results in relation to the semi-classical and 
quantum analyses which have been described above. The results confirm the emergence 
of a quantum phase boundary in the limit of zero atomic imbalance. 

2 The model 

We consider a general three-mode Hamiltonian describing a heteronuclear molecular Bose- 
Einstein condensate with two distinct species of atoms, labelled by a and b, which can com- 
bine to produce a molecule labelled by c. We introduce canonical creation and annihilation 
operators {a, b, c, a\ , c^} satisfying the usual commutation relations [a, a^} = I etc., 
which represent the three degrees of freedom in the model. The Hamiltonian reads [19] 

H = U aa N 2 a +U bh N 2 b +U CC N 2 C +U ab N a N b + U ac N a N c + U bc N b N c 

+ n a N a + ^ b N b + fi c N c + tt(a f b j c + c j ba). (1) 

The parameters Uy describe S-wave scattering, /ij are external potentials and O is the 
amplitude for interconversion of atoms and molecules. We remark that in the limit U aa = 
U bb = U cc = U ab = U ac = U bc = 0, equation (JTJ) is the Hamiltonian studied in [20,21] in 
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the context of quantum optics. In the latter stages of the manuscript we will study this 
limiting case in some detail. 

The Hamiltonian acts on the Fock space spanned by the (unnormalised) vectors 

\n a ,n b ;n c ) = (ay«(tf) n >>(c^\0) (2) 

where |0) is the Fock vacuum. We then have 

N a \n a ; n b ; n c ) = n a \n a ; n b ; n c ) 

etc., where N a = a^a, N b = tfb and N c = c^c. The Hamiltonian commutes with J = 
N a — N b and the total atom number N = N a + N b + 2N C . We refer to J as the atomic 
imbalance and introduce k = J/N, k G [—1,1] as the fractional atomic imbalance. As 
there are three degrees of freedom and three conserved operators, the system is integrable. 
This fact will allow us to analyse the model in some depth. Below we begin with a semi- 
classical analogue of the model, and determine the fixed points of the system. 



3 Semi-classical analysis 

Let Nj, <f>j, j = a, b, c be quantum variables satisfying the canonical relations 

[(f) j, <j> k ] = [Nj, N k ] = 0, [Nj, <f> k \ = i5 jk I. 

We make a change of variables from the operators {j, j' \ j = a, b, c} to a number-phase 
representation via 

j = eyv(i ( Pj)\fNj 3 = a, b, c 

such that the canonical commutation relations are preserved. We now make a further 
change of variables 

z = ^(N a + N b -2N c ), 
N 

= ~^{<Pa + <Pb - 0c), 

such that z and are canonically conjugate variables; i.e. 

[z, 0] = il. 

For large iV we can now approximate the (rescaled) Hamiltonian by 

H = \z 2 + 2(a - \)z + A - 2a + (3 + y/2(l - z)(z + c + )(z + cTj cos (^j (3) 

with 

V2N /Ugg Ubb Ucc Uab Ugc U b c 

QV4 4 4 4 4 4 
\/2N fl + k TT l-k TT 1 TT l + k TT l~k TT 1, 

a = —^-\-^-U aa + ^-U b b+-U ab — U ac — U bc+ ^{Va + ~ He) 

J2N { 2 
P = ^f(l + A;) 2 ^ a + (l-A;) 2 ^ + (l-A; 2 )f/ afe + -((l + A;)/i a + (l-A;K 



-z)(z + c + )(z + c-) sin f^j 



where c± = 1 ± 2fc. Since iV and are conserved, we treat them as constant. 

We now regard @ as a classical Hamiltonian and investigate the fixed points of the 
system. The first step is to derive Hamilton's equations of motion yielding 

dz _dH 4 

dt d(j) 

= ■» ■ (l-»)(2» + 2)-(» + c t )(, + 0.„/# 

A & V2(l-j)(z + c + )(; + c_) V'V 

The fixed points of the system are determined by the condition 

dH dH 

Due to periodicity of the solutions, below we restrict to G [0, Nn/2). It is necessary to 
treat the cases of k ^ and k = separately, and without loss of generality we assume 
k > 0. 

3.1 Case I: k ^ 

Define the functions 

f( z ) = \ z + a -\ (5) 
= (^-l)(2^ + 2) + (z + c + )(^ + C-) 
2^2(1 - 2 )(^ + c + )(z + c_) 

Note that the domain of g(z) is z e [2/c — 1,1], and is divergent at 2 = 2k — 1 and 
z = 1. For fc 7^ 0, we then have the following classification of solutions for (j3J): 

• = 0, and z is a solution of 

f(z)=g(z) (7) 
which can admit one, two or three solutions. 

• = Nn/ '4, and z is a solution of 

/(*) = (8) 
which can admit one, two or three solutions. 

A graphical representation of possible types of solutions for 9 = is given in Fig. 
^ From the equations (0 (HJ) we can determine there are fixed point bifurcations for 
certain choices of the coupling parameters. These bifurcations allow us to divide the 
coupling parameter space into different regions. To construct this diagram, we observe 
that bifurcations occur when / is the tangent line to g±; i.e. for values of A, a such that 



A = ± d -f 

dz 



(9) 

f(zo) = ±g(za) (10) 



for some Zq. This requirement determines the boundaries in parameter space, which are 
depicted in Fig. |21 
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Figure 1: On the left, graphical solution of (J7J) with k = 0.8. Depending on the values 
of A and a, there may be one, two or three solutions. On the right, graphical solution 
of (fT^j) with k = 0. Depending on the values of A and a, there may be zero, one or two 
solutions. 



3.2 Case II: k = 

Next we consider the case k = for which the function g(z) has substantially different 
properties. Setting c + = c_ = 1 into (jBJ), we find that g(z) reduces to 

/ % 1 3,2- , . 

g( z ) = — ; (11) 

2v/2(r^) 



Here we observe that g(z) is divergent at z = 1, but finite at z = — 1. This property affects 
the types of solutions for (HJ). Specifically, we now have the following classifications of 
solutions for k = 0: 

• = 0, and z is a solution of 

f(z) = g(z) (12) 

which can admit zero, one or two solutions. 

• <p = Nit /4, and z is a solution of 

f(z) = -g(z) (13) 

which can admit zero, one or two solutions. 

• z = — 1 and is a solution of 

cos J =-2X + a (14) 



which can admit zero, one or two solutions. 
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Figure 2: Parameter space diagram identifying the different types of solutions for equation 
(j3J), for k = 0.002, 0.02, 0.2. In each case the diagram is divided into regions A (one 
solution for z when = and one solution when = Ntt/4 ), B (three solutions for z 
when = and one solution when = Ntt/4) and C (one solution for z when = 
and three solutions when = Ntt/ '4). The boundary separating the regions is given by 
solutions to the equations (j9H0|) . 

A graphical representation of possible types of solutions for 9 = is given in Fig. 
[0 Because g{— 1) is finite, for this case there can be either zero, one or two solutions. 
As in the k ^ case, we can determine the region boundaries in parameter space from 
equations (j9llOJ) . Moreover, because of the existence of solutions of the form given by 
([T4|) for k = 0, which do not have an analogue for k ^ 0, we see the appearance of new 
boundaries given by the conditions A = (a ± l)/2 for all values of a. The boundaries in 
parameter space are depicted in Fig. El 

To help visualise the classical dynamics, it is useful to plot the level curves of the 
Hamiltonian ©• Since the fixed point bifurcations change the topology of the level 
curves, qualitative differences can be observed between each of the regions. The results 
are shown respectively in Fig. 0]for k = 0.2 and Fig. |S]for k = 0, where for clarity we 
show 40/iV G [-2?r, 2tt]. 

Hereafter we will focus most attention on the case where A = 0, so the model has one 
effective coupling parameter, a. From Figs. 121 H3 it can be seen that for this submani- 
fold there are no bifurcations when the atomic imbalance is non-zero, with bifurcations 
occuring at a = ±1 when the atomic imbalance is zero. For the case when the atomic 
imbalance is non-zero, the global minimum of the classical Hamiltonian (j3J) occurs when 

= Ntt/4: and z is the unique solution of ©• In particular, for the solution z G [2k — 1, 1] 
dz/da is a continuous function of a. When the atomic imbalance is zero and a > 1, the 
global minimum of the classical Hamiltonian Q always occurs at the phase space bound- 
ary z = — 1 with arbitrary. At a = 1 a bifurcation occurs, and for a slightly less than 

1 two saddle points arise for z = — 1 with given by solution to (|T3jl and a new global 
minimum emerges corresponding to = Ntt/4 with z the unique solution of (j!3|) . In this 
case dz/da is discontinuous at a = 1. 
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Figure 3: Parameter space diagram identifying the different types of solutions for equation 
(0J) when k = 0. In region I there is no solution for z when = 0, and one solution for z 
when = Nit/ 4. In region II there are two solutions for z when = 0, and one solution 
for z when = Nn/4. In region III exists one solution for z when = 0, one solution for z 
when = Nn/4, and two solutions for when z = —1. In region IV there is one solution 
for z when = 0, and no solution for z when = Nn/4. In region V there is one solution 
for z when = 0, and two solutions for z when = Ntt/4. The boundary separating 
regions II and III is given by A = (a + 1)/2, while the equation A = (a — 1)/2 separates the 
regions III and IV. The boundary between regions I and II has been obtained numerically. 

In the following sections we will conduct an analysis of the quantum Hamiltonian 
(JTJ). In particular we will establish that the bifurcation occuring at (a, A) = (1,0) when 
the atomic imbalance is zero can be seen to influence the ground-state properties of the 
quantum system. In the context of the quantum system we will refer to the boundaries in 
Figs. EHHlas threshold couplings. We avoid using the terminology quantum phase transition 
as the analysis is conducted for finite particle number, not in the thermodynamic limit. 
The reason for not taking the thermodynamic limit is that the quantities A and a are 
dependent on N. Additionally, in the thermodynamic limit iV — > oo with k finite the 
semi-classical results predict qualitative differences between the cases k = and k ^ 0. 
However if N is odd then we cannot have k = 0, raising technical issues about whether the 
limit is convergent. Consequently we only consider the case of finite particle number. To 
deal with the subtleties of the finite size of the system we will formally define a quantum 
phase pre-transition in Sect. El 

4 Exact solution of the quantum Hamiltonian 

We now turn our attention to a quantum mechanical treatment of the model, to investigate 
the nature of the additional threshold couplings when the atomic imbalance is zero. First 
we derive an exact Bethe ansatz solution of the model, and then use this to map the 
spectrum of the Hamiltonian into the spectrum of a one-body Schrodinger operator. 
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Figure 4: Level curves of the Hamiltonian (J3J) for k = 0.2, where the dark regions indicate 
lower values than the light regions. Figures (a) and (d) correspond to region A while 
Figures (b), and (c) correspond to region B. The parameter values are: (a) A = 10, a = 
—5; (b) A = 10, a = 4; (c) A = 10, a = 12 and (d) A = 10, a = 16. In region A, there 
is a maximal point at = and minima at 40/iV = ±7r. Two additional fixed points, a 
saddle and a maximum, occur in region B at <fi = 0. 

4.1 Energy eigenvalues as roots of a polynomial equation 

We rewrite the system Hamiltonian in a compact form as 



where the operator U is a function of the number operators: 

U = U aa N 2 a + U bb N 2 b + U CC N 2 C + U ab N a N b + U ac N a N c + U bc N b N c 
+ fi a N a + [i b N b + fi c N c . 

Since the operators iV and k = J/N are conserved we fix these and without loss of 
generality consider cases where k > 0. This restricts the Hilbert space to a subspace of 
dimension (m + 1) spanned by the vectors 



H = U + n(a ] tfc + c f 6a) 



(15) 



l-j;m-j;j) 



(16) 



where we have defined 



_ N(l + k) 

2 ' 

Ml - k) 

m = 

2 

such that I + m = N. We then look for eigenstates of of the form 



ni 




(17) 



j=0 
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Figure 5: Level curves of the Hamiltonian © for k = 0, showing the typical behaviour 
for regions I, II, III and IV. The dark regions indicate lower values than the light regions. 
The parameter values are A = 1.0, a = —2.0 for region I, A = 2.0, a = 2.0 for region II, 
A = 0.5, a = 0.5 for region III and A = 0.5, a = 3.0 for region IV. In region I there are 
local minima for 40/iV = ±7r. Besides the minima at 40/iV = ±7r, two additional fixed 
points (a maximum and a saddle point) are apparent in region II occurring at <p = 0. 
In region III there are minima at 4<ft/N = ±7r and for = just one fixed point, a 
maximum. There are also saddle points for when z — — 1. In region IV just one fixed 
point, a maximum, occurs for = 0, which always has z < 1. In contrast the global 
minimum occurs for z — —1. 

Since the basis states (fTp^) are eigenstates of each of the number operators they are 
eigenstates of the operator U so we can define the quantities Uj through 

U\l - j;m- j;j) =U j \l-j;m-j;j). (18) 

The Hamiltonian acts on the general state (|T7j) as 

m— 1 

H\^) = J2( U jPj + + !)Pi+i + (l + l-j)(m + l- i)pj-x)) \l -j;m- 

3=1 

+ (U po + Mpi) \l] m; 0) + {U m p m + Qp m -i(l -m+l))\l- m; 0; m) (19) 

Requiring that ()17|) is an eigenstate of the Hamiltonian with energy eigenvalue E leads 
to the following recursion relations that must be satisfied by coefficients pj\ 

n Pl +U po = E Po , (20a) 
ft((j + l)p j+1 + (l + l-j)(m + l- j) Pj _ x ) + UjPj = E Pv (20b) 

U m p m + Q(l - m + l)p m -i = Ep m , (20c) 

where 1 < j < in — 1 in (j20b|) . As the normalisation of the state (|17|) can be chosen 
arbitrarily, we have the freedom to choose po = 1. The recursion relation (j20b|) then 
shows that pj is a polynomial in E of order j. The constraint ()20c|) is thus a polynomial 
in E of order (m + 1), whose roots are the energy eigenvalues of Since the number of 
roots is the same as the dimension of the subspace spanned by the vectors (HBJ), all energy 
eigenvalues are given by the roots of (|20cj) . 
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4.2 Bethe ansatz solution and mapping to a Schrodinger equa- 
tion 



With the above implicit form for the energy eigenvalues we are able to map the energy 
spectrum into that of a one-dimensional Schrodinger equation. We start by mapping the 
energy eigenstates to polynomial solutions of a particular second-order ordinary differ- 
ential equation (ODE) and then utilise a change of variables such that the differential 
equation takes the form of the Schrodinger equation. Each eigenstate of the system (fT7j) 
can be represented by an m th order polynomial with coefficients pj (j = 1,2, ..m.). For 
a particular energy, we can then construct an ODE for G(u) such that the polynomial 
coefficients must satisfy the recursion relations of (ffilj) . Below we outline the details of 
this construction. 

Consider a general second-order ODE eigenvalue problem satisfied by an m th polyno- 
mial function G(u): 

a{u)G" + b(u)G' + c{u)G = EG (21) 
First we write the polynomial G(u) with roots {m p }^1 1 in the factorised form 

m 

P =i 



such that 



p=l q^p 



U q ), 



m m m 



p=l q^p r^p 
r+q 



Evaluating (|2*Tj) at the root u q leads to the Bethe ansatz equations 

b{u q ) 



:V , g = l,2,...,m. (22) 

Hence, the roots of the polynomial must satisfy the system of coupled equations if 
G{u) is a solution to (|2*T|). 

We can map the solutions of ()21|) with eigenvalue E to solutions of a Schrodinger 
equation 

~ d y^ + V{x)iP(x) = EiP(x) (23) 
ax z 

with the same eigenvalues, by mapping the polynomial solution of (|2ip to a wavefunction 
of (J23J) via 

$(x) = e f{x) G(u(x)). 
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Substituting into the Schrodinger equation gives the following relations to be satisfied 

.(«) = -(§' (24a) 

2 (-) 

^) = V( X )-^l-(^-) 2 (24c) 



dx 2 \ dx r 

In view of the above discussions, we now formulate the Bethe ansatz solution for 
and the associated mapping to a Schrodinger equation. To simplify the notation, we 
define 

Uj = A(m - j)(m - j - 1) + B(m - j) + C 



where 



A = U aa + Utb + U cc + Uat - U ac - U bc 

B = (l + 2l- 2m)U aa + U bb + (1 - 2m)U cc + (1 + I - m)U ab 

+ (2m - / - l)U ac + (m - l)U bc + fi a + fi b - fi c 
C = (1 - m) 2 U aa + m(Z - mjC/oc + m 2 [/ cc + (m - /)/i a + m/i c . 

The polynomial defined as 

m 

G(u) = Y,P^ m - 3 , (25) 

with the pj satisfying (|2(Jal2(Jbl2(Jc|) . is a solution to the following differential equation 

{Au 2 + Qu)G" + (Bu + Q(l - m + 1 - m 2 ))G' + (Qmu + C)G = EG. (26) 
The roots of G{u) are solutions of the Bethe ansatz equations 
n(l - m + 1 - u 2 ) + Bu q ™ 2 

— ro I 4 ^ = E > g = l,2,...,m. 27 

Mg(fi + AM,) ^ M p - 

We can also derive an expression for the energy eigenvalues of the model in terms of the 
roots u q . Consider the leading order expansions 

m 

G(u) =u m -u m ~ 1 ^u q + ... 

q=l 



G\u) = mu m - 1 - (m - l)u m - 2 ^u q + ... 



9=1 



G"{u) = m(m - l)u rn - 2 - (m - l)(m - 2)M m ~ 3 + ... 

<r=i 
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We substitute these expressions into ([26]) and equate terms of order m to arrive at the 
following expression for the energy eigenvalues of the system 



m 



E = Am(m - 1) + Bm + C - n^u q (28) 

Next we determine the explicit form of the Schrodinger equation. Comparing (J2T?|) to 
(ED) gives 

a{u) = Au 2 + flu 

b(u) = (I - m + 1 - u 2 )n + Bu 

c(u) = muQ + C 

Using (|24a|24bl24cJ) we may perform the mapping to the Schrodinger equation by choosing 

du 



— = ±V-Au 2 - Qu 
ax 

Integrating this expression (with a convenient choice for the constant of integration) gives 

u = ^(cos(v / Ix) - 1) (29) 

We also find that 

df Q 2 ( B Q 2 \ 

-j- = — s sm(VAx) + y/A(l - m + 1) = - — csc(VAx) 

dx 4A§ V 2VA 2A 3 / 2 J 



\fA VI 2 B \ /—— 

— + 2Av- 2 + m coti) 



So the wavefunction 



m ( Q 

= exp (/(x)) JJ ( — (cos(v / Ix) - 1) - 1 
p=i ^ 

satisfies the Schrodinger equation with potential 

V{x) = mvn + C+Pf+(!f) 
dx z \dx J 

/ Q 2 , . Q A A (\ 3Q 2 B 

fl 4 . 2/A7 , VI 2 ( tt 2 B\ . r- , 

+ Sm (v ^ } + 2A I" + 2^ + 2A ) C0S{VJX) 

f3A .„ 1N2 ft 2 n 1N ft 4 Q 2 n (B Q 2 

+ {T + A{l ~ m + 1) + ^ l - m + 1) + ^-^ + B {2A + A 2 

x csc 2 (VAx) 

+ (u +B J (/ - m+2) - 2A(/ - m+l) -2^- 5 U + ^j 



(30) 



X 



cot(v^Ax) csc(v^4x). 
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The above potential is an example of a quasi-exactly solvable potential [22], whereby a 
finite number of eigenstates of the form (|3U|) can be constructed. The concept of mapping 
the spectrum of many-body systems into those of one-body Schrodinger equations has 
been discussed in detail in [23]. 



5 Analysis in the no scattering limit 

In this section we now conduct a deeper analysis of the Hamiltonian in the no scattering 
limit where Ujk = for all j, k — a,b, c. In this limit the model simplifies substantially, 
yet remains sufficiently non-trivial to enable us to gain an understanding of the quantum 
behaviour through the Schrodinger equation mapping, ground-state expectation values 
and quantum dynamics. Specifically the no scattering limit corresponds to the coupling 
A = in the semi-classical analysis of Sect. El With reference to Fig. El there are two 
threshold couplings in the case of zero atomic imbalance. One occurs at (a, A) = (1,0), 
signifying the bifurcation of the global minimum of the Hamiltonian, while the other 
occurs at (a, A) = (—1, 0), signifying the bifurcation of the global maximum. In contrast 
there are no bifurcations along the line A = in Fig. El We focus our attention to the 
coupling (a, A) = (1, 0) as the bifurcation of the fixed point in phase space is associated 
with the ground state of the quantum system. 



5.1 Schrodinger equation mapping 

For small values of A, we can take series expansions for the trigonometric functions in the 
potential V(x) and wavefunction *f>(x). Then taking the limit A — > (corresponding to 
A = for the analogous classical system 0) the Schrodinger potential becomes 

V(x) = C-|(iV+l)+(j 2 -^x- 2 

B 2 Q 2 ,„ _A o Q 2 B 4 Q 4 



+ (N + 2))x 2 + x 4 + x b (31) 

V 16 8 / 32 256 K J 

where we now parametrise the system in terms of the variable J = I — m and N = I + m. 
Now consider a simple subclass of the general Hamiltonian 

H = fiN c + VL((Jb ] c + Jba) (32) 

We have mapped the general model to a Schrodinger equation in the previous section. 
The above case (J3*2|) corresponds to A = 0, B = — /z and C = m\i in equation (|3*T|) . The 
energy eigenstates map to solutions of the Schrodinger equation with potential 

v(x) = 40+21 + (> _ I) 

+ l(^-2^(iV + 2))^-^ + ^« (33) 
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The associated wavef unction is given by 

*(*) = xW*> exp (^f + tf\ f[ {^f - u p ) (34) 

with energy eigenvalues 

(N-J)/2 

E = —Q u i ( 35 ) 

9=1 

where the {u q } are solutions to the Bethe ansatz equations 




2 4 6 8 10 12 1 2 3 4 5 6 7 8 9 



x x 

(a) (b) 

Figure 6: The potential V(x) given by (a) iV = 501 and atomic imbalance J = 1. 

The potential is bounded from below, with the inset showing V — > oo as x —>■ 0. Varying 
the coupling parameter a across the threshold value a = 1, it is apparent there is no 
bifurcation of the potential minimum, (b) iV = 500 and atomic imbalance J = 0. The 
potential is not bounded below, with the inset showing V — > — oo as x — > 0. As the 
coupling parameter a is varied across the threshold value a = 1, it is apparent there is a 
bifurcation with the formation of a local minimum and a local maximum for a < 1. 

The semi-classical analysis predicts a threshold coupling at a = 1 when the atomic 
imbalance is zero where for the case under consideration a = —fi/(Qy2N). When the 
atomic imbalance is non-zero there is no predicted threshold coupling. Fig. El (a) depicts 
the potential ()33|) for N = 501, J = 1 and various values of a close to the threshold value 
a — 1. It can be seen that the potential has a single minimum for all a. In contrast, 
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Fig. IH1 (b) shows (}3"3"|) for iV = 500 and J = 0. For this case the potential is not bounded 
from below and there is a bifurcation for a ~ 1. For the model (}3"2j) . the predictions 
of the semi-classical analysis conducted in Sect. 01 of a threshold coupling at a ~ 1 are 
consistent with qualitative differences of the associated Schodinger equation. 

Now we examine bifurcations of the critical points of the potential (|33|) in more detail. 
Consider the general class of potentials 

V(x) = Ax~ 2 + Bx 2 + Cx 4 + Vx 6 (37) 

where C, T> are assumed to be positive and no constraints are placed on A nor B. In 
particular we wish to determine when the condition 

dV _ d 2 V _ 
dx dx 2 

can be met. Since the potential is a symmetric function, we restrict to non- negative values 
of x. When A = it is straightforward to deduce that, for any values of C and V, (JHHj) 
holds at x = when B = 0. For non-zero values of A we find 

^ = -2Ax~ 3 + 2Bx + ACx 3 + QVx 5 (39) 
dx 

d 2 V 

-— = 6Ax~ 4 + IB + YlCx 2 + 30£>x 4 (40) 

cur 

We set both (|39|) and ()40j) to zero and take particular linear combinations to obtain the 
following relations: 

'S+l^ - («) 

*^_5*W = ^_^_ &6 = (42) 
8 cfe 2 8 v 7 

Note that eq. (p£T|) is independent of A, and has solutions 

2 -3C ± V9C 2 - 2ABV 

x = — • (43) 

12V v ; 



We take the positive square root in (|4*3^1 and impose B < 0, ensuring x 2 > 0. Treating i3 
as a small parameter such that 



yields 



l*| « % (44) 



B 

x 2 w (45) 



3C 

Next we substitute (|4*Hj) into (jl2|) and solve for B: 

B = 3(AC 2 ) 1/3 . (46) 
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Since B is negative such a solution only exists when A is negative. 
Matching the co-efficients between ()33|) and (pTTJ) gives 

A = J 2 - 1 - 
4 

B = 4(^-2^(^ + 2)) 



16 
/zfi 2 

32" 

n 4 



256 



When a ~ 1, or equivalently fi w — fiv2iV, we satisfy the requirements C, T> > and 
(jSJ). A bifurcation of the potential only occurs when the atomic imbalance J is zero, as 
A must be negative. Using (j^Sj) we then find satisfies 

3/i 2 / 3 fi 4 / 3 + /i 2 = 2ft 2 (iV + 2). 

From this expression we determine the leading quantum correction to the semi-classical 
result for the threshold coupling for (|32j) : 

H w -fi(2(JV + 2)) 1 / 2 + : y(2(AT + 2))- 1 / 6 . 
5.2 Ground- state expectation values and quantum dynamics 




-30 -28 -26 -24 -22 -20 -18 -16 -14 -12 



Figure 7: Ground-state expectation values of the (scaled) molecular number operator N c 
as a function of the coupling fi, for the Hamiltonian ([32)1 . Results shown correspond to 
f2 = 1 and both zero and non-zero atomic imbalance. The inset shows the first derivative 
of the expectation values with respect to the coupling fi. While there are quantitative 
differences there is no significant qualitative change between the case of zero and minimal 
non-zero atomic imbalance. 
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Next we examine the behaviour of the ground state of (|32|) as the threshold coupling 
a = 1 is crossed. From the semi-classical analysis we have found that the global minimum 
for a > 1 and J = occurs at z = —1 in phase space. In the Hilbert space of states this 
corresponds to |0; 0; N/2). It is then appropriate to compute the gound-state expectation 
value (N c ) for the quantum system as the coupling is varied. Results are shown in Fig. 
|7| In general agreement with the semi-classical result, it can be seen that the expectation 
value (N c )/N is close to unity when — > Q\r2N for the case of zero atomic imbalance. 
When —jj, < Qy2N the expectation value decreases. The figure also shows that when 
the imbalance J = 1 the results are qualitatively similar. However from the predictions 
of both the semi-classical analysis of Sect. Eland the the associated one-body Schrodinger 
potential of Sect. 15.11 we do not obtain any prediction about the change in the ground- 
state properties when the imbalance is non-zero. Further increase in the value of J (not 
shown) does not indicate any dramatic change in the qualitative features of (2N C ) /N 
as a function of J. As mentioned previously, because we are studying a finite system 
changes in the ground-state properties are smooth as J is varied. If we instead look at 
the quantum dynamics as the threshold coupling is crossed, qualitative differences are 
more apparent. 



N=500 (1=1 J=0 N=500 (1=1 J=10 

H = H b =0 u.= li>=IO,0,250> |l,= HfO 11=0 Ii>=ll0,0,245> 




Figure 8: Time evolution of the expectation value of z for the Hamiltonian (}3*2*)) with 
N = 500. The cases shown are, from top to bottom, a = 0.9, 0.95, 1, 1.05, 1.1. (a) 
J = and initial state |0; 0; 250). The oscillations are largely irregular with significantly 
decreasing amplitude as the point at a = 1 is crossed. This point corresponds to the 
boundary at (a, A) = (1, 0) between regions III and IV as shown in Fig. (b) J = 10 
with initial state 1 10; 0; 245). The oscillations display collapse and revival behaviour with 
smoothly decreasing amplitude as the point at a = 1 is crossed, indicative of the fact 
there is no boundary at (a, A) = (1, 0) in Fig. El 
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In general the time evolution of any state is given by \^{t)) = U(t)\<p), where U(t) 

m 

is the temporal evolution operator U{t) = \ ex p( — iEjt) , \j) is an eigenstate with 

j=o 

energy Ej and \<p) represents the initial state with N = N a + Nf, + 2N C . We adopt the 
method of directly diagonalising the Hamiltonian ()32j) as done in [18], and compute the 
expectation value of z(t) through 

{z(t)) = ±(*{t)\N a + N b -2N e \*(t)). 

For a fixed atomic imbalance J we will use the initial state configuration | J; 0; (N — J)/2). 
When J = this state correspond to z = —1 in phase space, which is a fixed point when 
a > 1. We thus expect that in this case (z(t)) will not vary significantly in time (i.e. 
the system is localised). On the other hand when J ^ the state | J; 0; (N — J)/ 2) does 
not correspond to a fixed point. We therefore compare the two cases of the quantum 
dynamics, with atomic imbalance J = and J ^ 0, as the value a = 1 is crossed. We fix 
the parameter Q = 1 and use /i as the variable coupling parameter. 

Results of the expectation value for z are shown in Fig. [H] for the cases of zero and 
non-zero atomic imbalance. The qualitative difference are quite apparent. In the case 
of zero atomic imbalance (k — 0), Fig. |H1 (a), we find that for a < 1 there are irregular 
oscillations in z. By comparison the dynamics in Fig. |H1 (b) for non-zero imbalance 
(k = 0.02) show a collapse and revival of oscillations. As the coupling parameter a is 
increased across the threshold value at a = 1, the transition to localised oscillations is 
much sharper in case (a) compared to case (b). Note in particular the vertical scales in 
(a) and (b) are not the same. We remark that the nature of the dynamics for a > 1 
does change smoothly from localisation to delocalisation over the intermediate values 
< k < 0.02 (not shown). Taking the thermodynamic limit iV — ► oo does not aid in the 
analysis. For the Hamiltonian (|32|) the condition for localisation of oscillations for zero 
atomic imbalance, a > 1, is equivalent to —fi > il\ / 2N. Hence for fixed —fi > and 
Q > this condition imposes an upper bound on N for which localisation occurs. 



6 Wavefunction overlaps 

In order to gain a better insight into the effect of the threshold couplings for the quantum 
system, in our final analysis we adopt the method of wavefunction overlaps [12,13]. If a 
system admits a quantum phase transition, then two states belonging to different phases 
of the same system are distinguishable. If states are distinguishable they must be orthog- 
onal [24] and consequently the wave function overlaps vanish. For systems which exhibit a 
quantum phase transition in the thermodynamic limit, the wavefunction overlaps between 
states in different phases go to zero in this limit. The occurrence of a minimum in the 
incremental ground-state wavefunction overlap in a finite system is then a precursor for a 
quantum phase transition in the thermodynamic limit. Thus for finite systems we iden- 
tify quantum phase pre-transitions at couplings for which the incremental wavefunction 
overlap is (locally) minimal. 
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Figure 9: Ground-state wavefunction overlaps for the Hamiltonian (|32j) with N = 1000 
and various values of A. The value of the local minimum at a ~ 1 is a decreasing function 
of A, asymptotically approaching zero. 

We now formally define a quantum phase pre-transition in terms of ground-state wave- 
function overlaps. Let H(5) denote a generic Hamiltonian depending on a coupling pa- 
rameter 5. Assuming the ground state of the system is non-degenerate, let 1^(5)) denote 
the unique normalised ground state. For fixed small A we define the function W&(5) by 

Wa{6) = \(*(6(1-A))\*{6{1 + A)))\ 

which is symmetric in A, bounded between and 1, and satisfies Wo (5) = 1- Generically, 
Wa(8) is a decreasing function of A. Fig. El shows the behaviour of the wavefunction 
overlaps for the Hamiltonian ()32)1 with N = 1000, and different values of A. It is clear 
that there is a distinct dip in the quantity Wa{o) near the threshold coupling a = 1. The 
choice for A affects the magnitude of the minimum, which can be made arbitrarily small. 
In Fig. El when A = 0.05, representing a coupling change of about 5% on either side of 
the threshold coupling, the ground states are essentially orthogonal. However the value 
of a at which the minimum occurs is largely independent of A. For a given A we say 
that there is a quantum phase pre-transition at 5 C if W&(5), treated as a single- variable 
function of 5, is locally minimal at 5 C . 

We have computed the wavefunction overlaps for several cases with both zero and non- 
zero atomic imbalance. Fig. (a) shows the behaviour of Wa(qz) with A = 0, A = 0.01 
and varying N for both k = and k = 0.02. It is clear the minimum value of Wa(q), 
which determines the quantum phase pre-transition, is at a ~ 1. The distinction between 
the predicted threshold coupling and the observed pre-transition coupling is that the pre- 
transition coupling also occurs for k ^ 0, although for fixed N the value of the minimum 
is lower for k = compared to k — 0.02. In all instances the value of the mimimum 
decreases with increasing N. Fig. El (b) shows similar results for fixed iV = 1000 and 
varying A. In this latter case we see that the occurences of the mimima, determining the 
pre-transition couplings, fit well with the predicted boundary of threshold couplings given 
by A = (a — l)/2 (cf. Fig. EJ). However we again observe pre-transitions couplings for 
k = 0.02 which were not predicted in Sect. H3 
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Figure 10: (a) Ground-state wavefunction overlaps of the Hamiltonian (J32*j) . for N = 
500, 1000, 1500. The solid lines correspond to cases when the atomic imbalance is zero, 
while the dashed lines illustrate the behaviour for the fractional imbalance k = 0.02. Two 
general properties that can be observed at the pre-transition coupling a ~ 1 are (i) the 
minimum value decreases with increasing N\ (ii) for fixed N, the value of the minimum 
is lower for k = compared to k ^ 0. (b) Ground-state wavefunction overlaps of the 
Hamiltonian ()32|) for iV = 1000 and different values of A. The solid lines correspond to 
cases when the atomic imbalance is zero, while the dashed lines illustrate the behaviour 
for the fractional imbalance k = 0.02. The locations of the minima fit the line of threshold 
couplings given by A = (a — l)/2 as predicted by the semi-classical analysis. 

In the above cases the minimum of Wa(oc) is substantially more pronounced for zero 
imbalance compared to non-zero imbalance. We can see that although quantum phase pre- 
transitions as defined exist for k ^ 0, the distiguishability of the two phases is more reliable 
in the limit as k goes to zero. In the analyses of Sect. Eland Sect. 15. II qualitative differences 
are only found precisely when k = 0. We interpret these results as the emergence of 
quantum phase boundaries at k — 0. 



7 Discussion 

We have studied the ground-state phases of a three-mode model describing a heteronu- 
clear molecular Bose-Einstein condensate, through a variety of techniques. Using a semi- 
classical analysis we were able to determine threshold couplings associated with fixed 
point bifurcations in phase space. We then derived the exact Bethe ansatz solution for 
the system and discussed how the spectrum of the Hamiltonian maps into that of an asso- 
ciated one-body Schrodinger equation. It was shown that in the particular subcase of no 
scattering interactions the threshold coupling for the global minimum in phase space for 
zero atomic imbalance was consistent with the existence of a bifurcation of the potential 
of the Schrodinger equation. For the non-zero imbalance case where the semi-classical 
results do not predict any threshold coupling, it was found that there was no bifurcation 
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of the Schrodinger potential. These results suggested that the ground-state properties of 
the model are sensitive to whether the atomic imbalance was zero or non-zero. However 
due to the finite nature of the system, such a transition is not associated with a discon- 
tinuity that is defined in the thermodynamic limit. We then introduced the notion of a 
quantum phase pre-transition for finite systems, defined in terms of ground-state wave- 
function overlaps. Applying this idea to the model under consideration we argued that a 
quantum phase boundary emerges in the limit as the atomic imbalance goes to zero. 

For future work it would be useful to investigate the ground-state entanglement prop- 
erties of the Hamiltonian. Not only would this be of interest in the study of the behaviour 
of the entanglement at the quantum phase boundary as k approaches zero, but the model 
is also a simple example of a strictly tripartite system which may offer insights into the 
role of three-way entanglement in the description of quantum phases. 
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